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CORONAL HEATING. WEAK MHD TURBULENCE AND SCALING LAWS 
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ABSTRACT 

Long-time high-resolution simulations of the dynamics of a coronal loop in cartesian geometry are 
carried out, within the framework of reduced magnetohydrodynamics (RMHD), to understand coronal 
heating driven by motion of field lines anchored in the photosphere. We unambiguously identify MHD 
anisotropic turbulence as the physical mechanism responsible for the transport of energy from the large 
scales, where energy is injected by photospheric motions, to the small scales, where it is dissipated. As 
the loop parameters vary different regimes of turbulence develop: strong turbulence is found for weak 
axial magnetic fields and long loops, leading to Kolmogorov-like spectra in the perpendicular direction, 
while weaker and weaker regimes (steeper spectral slopes of total energy) are found for strong axial 
magnetic fields and short loops. As a consequence we predict that the scaling of the heating rate 
with axial magnetic field intensity Bq, which depends on the spectral index of total energy for given 

loop parameters, must vary from B^ 2 for weak fields to Bq for strong fields at a given aspect ratio. 
The predicted heating rate is within the lower range of observed active region and quiet Sun coronal 
energy losses. 

Subject headings: Sun: corona — Sun: magnetic fields — turbulence 



1. INTRODUCTION 

In this letter we solve, within the framework of RMHD 
in cartesian geometry, th e Parker field-line tangling 
(coronal heating) problem (|Parkerlll972l 119881 ). We do 
this via long simulations at high resolutions, introduc- 
ing hyper-resistivity models to attain extremely large 
Reynolds numbers. We show how small scales form and 
how the coronal heating rate depends on the loop and 
photospheric driving parameters, and derive simple for- 
mulae which may be used in the coronal heating context 
for other stars. 

Over the years a number of numerical experiments have 
been carried out to investigate coronal heating, with par- 
ticular emphasis on exploring how photospheric field line 
t angling leads to cur rent sheet formation. 

iMikic et all (fl989l) and iHendrix fc Van Hovenl (fl996l ) 
first carried out simulations of a loop driven by pho- 
tospheric motions using a cartesian approximation (a 
straightened out loop bounded at each end by the pho- 
tosphere) imposing a time-dependent alternate direction 
flow pattern at the boundaries. A complex coronal mag- 
netic field results from the photospheric field line ran- 
dom walk, and though the field does not, strictly speak- 
ing, evolve through a sequence of static force-free equi- 
librium states (the original Parker hypothesis), mag- 
netic energy nonetheless tends to dominate kinetic en- 
ergy in the system. In this limit the field is struc- 
tured by current sheets elongated along the axial direc- 
tion, separating quasi-2D fl ux tubes which con s tantly 
move around and interact. lLongcope fc Sudani (|1994T ) 
focused on the current sheet formation process within 
the RMHD approximation, also used in the simulations 
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by iDmitruk fc G6me3 (|1999f ). The results from these 
studies agreed qualitatively among themselves, in that all 
simulations display the development of field aligned cur- 
rent sheets. However, estimates of the dissipated power 
and its scaling characteristics differed largely, depending 
on the way in which extrapolations from low to large 
values of the plasma conductivity of the properties such 
as inertial range power law indices were carried out. 2D 
numerical sim ulations of incompressible MHD with mag- 
netic forcing ([Einaudi et al.lll99dlGeorgoulis et ailll998t 
IDmitruk et al.lll998HEinaudi fc Vellilll999h showed that 
turbulent current sheets dissipation is distributed inter- 
mittently, and that the statistics of dissipation events, 
in terms of total energy, peak energy and event duration 
displays power laws not unlike the distribution of ob- 
served emission events in optical, ultraviolet and x-ray 
wavelengths of the quiet solar corona. 

More recently full 3D sections of the solar corona 
with a realistic geometry have been simulated by 
Gudi ksen fc Nordlundl (|2005f ). While this approach has 
advantages when investigating the coronal loop dynam- 
ics within its neighboring coronal region, modeling nu- 
merically a larger part of the solar corona drastically 
reduces the number of points occupied by the coronal 
loops. Thus, these simulations have not been able to 
shed further light on the physical mechanism responsible 
for the coronal heating. 

In § [H we introduce the coronal loop model and the 
simulations we have carried out; in § [3] we describe our 
numerical results, and in §[4] we give simple scaling argu- 
ments to understand the magnetic energy spectral slopes. 
This will lead to a quantitative asymptotic estimate of 
the coronal loop heating rate, and of its scaling with the 
axial magnetic field, photospheric velocity amplitude and 
coronal loop length. 

2. THE MODEL 

A coronal loop is a closed magnetic structure threaded 
by a strong axial field, with the footpoints rooted in 
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Fig. 1. — : High-resolution simulation with vj^/u p h = 
200, 512x512x200 grid points and H x = 800. Magnetic 
(Em) and kinetic (Ek) energies as a function of time 
(T4 = L/v_a is the axial Alfvenic crossing time). 

the photosphere. This makes it a strongly anisotropic 
system, as measured by the relative magnitude of the 
Alfven velocity va ~ 1000 kms -1 compared to the typ- 
ical photosphcric velocity u p h ~ 1 kms -1 . This means 
that the relative amplitude of the Alfven waves that are 
launched into the corona is very small. The loop dynam- 
ics may be studied in a simplified geometry, neglecting 
any curvature effect, as a "straightened out" cartesian 
box, with an orthogonal square cross section of size £±, 
and an axial length L embedded in an axial homogeneous 
uniform magnetic field Bo = i?o e z- This system may 
be described by the reduced MHD (RMHD ) equations 
(|Kadomtsev fc Pogutsd Il974t IStrausd Il976h : introduc- 
ing the velocity and magnetic field potentials (p and tp, 
u± = V x (<pe z ), b± — V x (tpe z ), and vorticity and 

current, uj = — V^i^s, j = — V^V-" the non-dimensioned 
RMHD system is given by 



dip dip 
duj dj 



f-lV n+1 



(1) 



(-1) 



n+1 



-V 2 ^lj. (2) 



As characteristic quantities we use the perpendicular 
length of the computational box £j_, the typical photo- 
spheric velocity u p h, and the related crossing time t± = 
£_i_/u p h- The equations have been rendered dimensionless 
using velocity units for the magnetic field (the density 
in the loops p is taken to be constant) and normaliz- 
ing by Uph. Then the non-dimensioned Alfven speed va 
in eqs. U])-© is given by the ratio va/u p u between the 
dimensional velocities. The Poisson bracket of two func- 
tions g and h is defined as [g,h] = d x gd y h — d y gd x h, 
where x,y are transverse coordinates across the loop 
while z is the axial coordinate along the loop. A simpli- 
fied diffusion model is assumed and lZ n is the Reynolds 
number, with n the hyperdiffusion index (dissipativity): 
for n = 1 ordinary diffusion is recovered. 

The computational box spans < x, y < 1 and 
< z < L, with L = 10, corresponding to an aspect 
ratio equal to 10. As boundary conditions at the photo- 



Fig. 2. — : Same simulation of Figure [T] The integrated 
Poynting flux S dynamically balances the Ohmic (J) and 
viscous (f2) dissipation. Inset shows a magnification of 
total dissipation and S for 200 < </r_4 < 300. 

spheric surfaces (z = 0, L) we impose a velocity pattern 
intended to mimic photosphcric motions, made up of two 
independent large spatial scale projected convection cell 
flow patterns. The wave number values k excited are all 
those in the range 3 < k < 4, and the average injection 
wavenumber is ki n ^3.4. 

3. RESULTS 

Plots of the rms magnetic and kinetic energies as a 
function of time, together with the dissipation due to cur- 
rents, vorticity, as well as the integrated Poynting flux, 
are shown in Figures Q] and [2J As a result of the photo- 
spheric forcing, energy in the magnetic field first grows 
with time, until it dominates over the kinetic energy by a 
large factor, before oscillating, chaotically, around a sta- 
tionary state. Fluctuating magnetic energy Em is ~ 35 
times bigger than kinetic energy Ek- 

The same generic features are seen in the rms current 
and vorticity dissipation, where however the time depen- 
dence of the signal is more strongly oscillating. The 
ohmic dissipation rate J is ~ 6.5 times viscous dissi- 
pation Q. The Poynting flux, on average, follows the 
current dissipation (there is no accumulation of energy 
in the box), however a detailed examination shows that 
the dissipation time-series tends to lag the Poynting flux, 
with notable de-correlations around significant dissipa- 
tion peaks. The spatial configuration of the currents 
which corresponds to a snapshot at a given time is dis- 
played in Figure [5] The currents collapse into warped, 
torn sheets which extend almost completely along the 
loop. The current peaks are embedded within the 2D 
sheet-like structures, corresponding to an anisotropic 
structure for the turbulence, in agreement with previous 
results. 

A dimensional analysis of eqs. (U)-© shows that the 
only free nondimcnsional quantity is / = £\_vaI Lu p h- 
We fix L/£± = 10 and vary the ratio of the Alfven 
speed to photospheric convection speed Vj\/u p h- Both 
runs with standard second order dissipation (n — 1) as 
well as hyperdiffusion (n = 4) have been carried out to 
obtain extended inertial ranges in the resulting spectra. 
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Fig. 3. — : Time-averaged total energy spectra for sim- 
ulations with v A /u ph = 50, 200, 400, 1000. Hyperdif- 
fusion (n = 4) has been used with ft 4 = 3 • 10 20 , 10 20 , 
10 19 , 10 19 respectively, and a grid with 512 x 512 x 200 
points. 

The power spectrum of total energy in the simula- 
tion box, once a statistically stationary state has been 
achieved, depends strongly on the ratio v A I u v h. Thi s was 
first found in simulations by iDmitruk et aLl (|2003[ ), de- 
voted to understanding how anisotropic regimes of MHD 
turbulence depend on boundary driving strength, with 
whom our numerical work is in broad agreement. 

The total energy spectrum, for values of v A /u p h = 
50, 200, 400, 1000 is shown in Figure H together with 
fits to the inertial range power law. As v A /u p h increases, 
the spectrum steepens visibly (note that the hump at the 
high wave- vector values for the runs with large v A /u p h is 
a feature, the bottleneck effect, which is well known and 
documented in spectral simulati ons of turbulenc e with 
the hyperdiffusion used here, e.g. lFalkovichl (|19 94)). with 
the slopes ranging from —2 to almost —3. At the same 
time while total energy increases, the ratio of the mean 
magnetic field over the axial Alfven velocity decreases, 
in good accordance with the theory. This steepening, 
which may be interpreted both as the effect of inertial 
line-tying of the coronal magnetic field and the progres- 
sive weakening of non-linear interactions as the magnetic 
field is increased, has a strong and direct bearing on the 
coronal heating scaling laws. 

4. DISCUSSION 

A characteristic of anisotropic MHD turbulence is that 
the cascade takes place mainly in the plane orthog- 
onal to the DC magnetic guide field (Shcbalin et al. 
1983). Consider then the anisot ropic version of the 
Iroshnikov-Kraichnan (IK) theo ry ([Sridhar fc Goldreichl 
11994 iGoldreichfe Sridharlll997f) . Dimensionally the en- 
ergy cascade rate may be written as p8z\ /T x , where 
8z\ is the rms value of the Elsasser fields z^ = tt^ ± 
at the perpendicular scale A, where because the system 
is magnetically dominated 5z^ ~ Sz^ . p is the average 



density and T\ is the energy transfer time at the scale A, 
which is greater than the ed dy turnover time t\ ~ \/5z\ 
because of the Alfven effect (Iroshnikov 1964; Kraichnan 
fl965l) . 

In the classical IK case, T\ ~ ta{t\/ta) . This cor- 
responds to the fact that wave-packets interact over an 
Alfven crossing time (with t\ > t a ), and the collisions 
follow a standard random walk in energy exchange. In 
terms of the number of collisions N\ that a wave packet 
must suffer for the perturbation to build up to order 
unity, for IK N\ ~ (t x /t a ) 2 . 

More generally, however, as the Alfven speed is in- 
creased the interaction time becomes smaller, so that 
turbulence becomes weaker and the number of collisions 
required for efficient energy transfer scales as 



with 



a > 2, 



(3) 



where a is the scaling index (note that a — 1 corresponds 
to standard hydrodynamic turbulence), so that 
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Integrating over the whole volume, the energy transfer 
rate becomes 



£±L- 



L 



L 



(5) 



Considering the injection scale A ~ £±, eq. © becomes 



e ~ l\L ■ p 



Sz] 



p£\L c 



(6) 



On the other hand the energy injection rate is given 
by the Poynting flux integrated across the photosphcric 
boundaries: ej n = pv A JdaUph ■ b±. Considering that 
this integral is dominated by energy at the large scales, 
due to the characteristics of the forcing function, we can 
approximate it with 



p£ 2 j_v A u ph Szi ± , 



(7) 



where the large scale component of the magnetic field can 
be replaced with Szi ± because the system is magnetically 
dominated. 

The last two equations show that the system is self- 
organized because both e and e,-„ depend on Sz£ ± , the 
rms values of the fields z at the scale £±: the inter- 
nal dynamics depends on the injection of energy and 
the injection of energy itself depends on the internal 
dynamics via the boundary forcing. Another aspect of 
self-organization results from our simulations: the per- 
pendicular magnetic field develops few spatial structures 
along the axial direction z, and in the nonlinear stage its 
topology substantially departs from the mapping of the 
boundary velocity pattern which characterizes its evolu- 
tion during the linear stage. T hese and other features will 
be discussed more in depth in lRappazzo et all (|2007f) . 

In a stationary cascade the injection rate ((7]) is equal 
to the transport rate ([6]). Equating the two yields for the 
amplitude at the scale £±: 



u p h V Luph ) 



(8) 
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Fig. 4. — : The solid line shows the exponent a/(a + l) as 
a function of a. Symbols show values of a corresponding 
to different values of v^/u p h, at fixed L/£± — 10. 

Substituting this value in ([6]) or (J7j) we obtain for the 
energy flux 



r ,2 2 ( Z-LVA X 



(9) 



where i>a = Bo/y/Aitp. This is also the dissipation rate, 
and hence the coronal heating scaling. The crucial pa- 
rameter here is / = £±v^/ Lu p h because the scaling in- 
dex a ([3]), upon which the strength of the stationary 
turbulent regime depends, must be a function of / it- 
self. The relative amplitude of the turbulence 8z^ ± /v^, 
is a function of /, and as / increases the effect of line- 
tying becomes stronger, decreasing the strength of tur- 
bulent interactions (wave-packet collision efficiency be- 
comes sub-diffusive) so that a increases above 2. The 
ratio 5z% jvj^ can also be interpreted as the rms value of 
the Parker angle Op, and is given by 



(10) 



This is actually an estimate of the average inclination of 
the magnetic field lines, while the rms value of the shear 
angle between neighboring field lines is at least twice that 
given by eq. (fT0|) , not considering that close to a current 
sheet an enhancement of the orthogonal magnetic field is 
observed (which leads to a higher value for the angle). 

Numerical simulations determine the remaining un- 
known nondimensional dependence of the scaling index 
a on /. The power law slopes of the total energy spectra 
shown in Figure [3] are used to determine a. Identify- 
ing, as usual, the eddy energy with the band-integrated 



Fourier spectrum Sz 
eq. ([5]) we obtain 



k± Ek ± , where k± ~ £±/X, from 



E kj _ (x k ± 



x + 2 



where for a = 1 the —5/3 slope for the "anisotropic 
Kolmogorov" spectrum is recovered, and for a = 2 the 
—2 slope for the anisotropic IK case. At higher values of 
a correspond steeper spectral slopes up to the asymptotic 
value of —3. 

In Figure |4] we plot the values of a determined in 
this way, together with the resulting power dependence 
a /(a + 1) of the amplitude ([8]) and of the energy flux ([9]) 
on the parameter /. The other power dependences are 
easily obtained from this last one, e.g. for the energy 
flux ([9]) the power of the axial Alfven speed v_a is given 
by 1 + a I (a + 1), so that in terms of the magnetic field 

Bo it scales as B^ 2 for weak fields and/or long loops, to 
Bq for strong fields and short loops. 

Dividing eq. © by the surface £\ we obtain the en- 
ergy flux per unit area F — e* j£\. Taking for ex- 
ample a coronal loop 40, 000 km long, with a num- 
ber density of 10 10 cm" 3 , v_\ — 2,000 kms -1 and 
u p h = 1 kms" 1 , (for these parameters we can estimate 
a value of a /(a + 1) ~ 0.95), which models an active 
region loop, we obtain F ~ 5 • 10 6 erg cm" 2 s" 1 and a 
Parker angle flTUJl < Op >~ 4°. On the other hand, 
for a coronal loop typical of a quiet Sun region, with a 
length of 100,000 km, a number density of 10 10 cm -3 , 
va = 500 kms" 1 and u p h = 1 kms" 1 , (for these param- 
eters we can estimate a value of a/(a + 1) ~ 0.7) we 
obtain F ~ 7 • 10 4 erg cm" 2 s" 1 and < © P >~ 0.9°. 

In summary, with this paper we have shown how coro- 
nal heating rates in the Parker scenario scale with coronal 
loop and photosphcric driving parameters, demonstrat- 
ing that field line tangling can supply the coronal heat- 
ing energy requirement. We also predict that there is 
no universal scaling with axial magnetic field intensity, a 
feature which can be tested by observing weak field re- 
gions on the Sun, or the atmospheres of other stars with 
differing levels of magnetic activity. 



(11) 



M.V. thanks W.H. Matthaeus for useful discussions. 
R.B.D. is supported by NASA SPTP. 



REFERENCES 



Dmitruk, P., Gomez, D. O., & DeLuca, D. D. 1998, ApJ, 505, 974 

Dmitruk, P., & Gomez, D. O. 1999, ApJ, 527, L63 

Dmitruk, P., Gomez, D. O., & Matthaeus, W. H. 2003, Phys. 

Plasmas, 10, 3584 
Einaudi, G., Velli, M., Politano, H., & Pouquet, A. 1996, ApJ, 457, 

L113 

Einaudi, G., & Velli, M. 1999, Phys. Plasmas, 6, 4146 
Falkovich, G. 1994, Phys. Fluids, 6, 1411 

Georgoulis, M. K., Velli, M., & Einaudi, G., 1998, ApJ, 497, 957 
Goldreich, P., & Sridhar, S. 1997, ApJ, 485, 680 
Gudiksen, B. V., & Nordlund, A. 2005, ApJ, 618, 1020 
Hendrix, D. L., & Van Hoven, G. 1996, ApJ, 467, 887 
Iroshnikov, P. S. 1964, Sov. Astron., 7, 566 

Kadomtsev, B. B., & Pogutse, O. P. 1974, Sov. J. Plasma Phys., 
1, 389 



Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385 

Longcope, D. W., & Sudan, R. N. 1994, ApJ, 437, 491 

Mikic, Z., Schnack, D. D., & Van Hoven, G. 1989, ApJ, 338, 1148 

Parker, E. N. 1972, ApJ, 174, 499 

Parker, E. N. 1988, ApJ, 330, 474 

Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B., ApJ, 
in preparation 

Shebalin, J. V., Matthaeus, W. H., & Montgomery, D. 1983, J. 

Plasma Phys., 29, 525 
Sridhar, S. & Goldreich, P. 1994, ApJ, 432, 612 
Strauss, H. R. 1976, Phys. Fluids, 19, 134 



Coronal Heating 



5 




Fig. 5. — : Top: side view of two isosurfaces of the squared current at a selected time for a numerical simulation with 
VA/uph = 200, 512x512x200 grid points and a Reynolds number IZi — 800. The isosurface at the value j 2 — 2.8 • 10 5 is 
represented in partially transparent yellow, while red displays the isosurface with j 2 — 8 • 10 5 , well below the value of 
the maximum of the squared current that at this time is j 2 — 8.4 • 10 6 . N.B.: The red isosurface is always nested inside 
the yellow one, and appears pink in the figure. The computational box has been rescaled for an improved viewing, 
but the aspect ratio of the box is 10, i.e. the axial length of the box is ten times bigger than its orthogonal length. 
Bottom: top view of the same two isosurfaces using the same color display. The isosurfaces are extended along the 
axial direction, and the corresponding filling factor is small. 



